function generateModelEulerError_file(allModelParams,x,y,xp,yp,strucOfArrays,positionFile,nameOfFunction)
% This function generates an M function which computes the Euler equation errors 
% for the model using either standard or extended perturbation.
% The results are saved in the file nameOfFunction, for instance 
% nameOfFunction = 'EulerEqError.m';

% We start deleting the old version of the file - if it exists
saveFile = [positionFile,'\',nameOfFunction];
if exist(saveFile,'file') > 0
    delete(saveFile)
end

text = ['function [ResEuler,yGrid,xGrid] = ',nameOfFunction(1:end-2),'(model,extendedPerOn,orderApp,setup)'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Some indices ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'ny      = model.ny;        %Number of output variables';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx      = model.nx;        %Number of state variables ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Unfold params';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(allModelParams)
    text = [allModelParams{i}, '= model.params.',allModelParams{i},';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end 
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The raw grid for the states ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if setup.xPointsDim == 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    xGrid = transpose(setup.xSim); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'else ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    dx       = (setup.upperx-setup.lowerx)/(setup.xPointsDim-1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    xGridtmp = zeros(nx,setup.xPointsDim); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    for j=1:setup.xPointsDim ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        xGridtmp(:,j) = setup.lowerx(:,1)+(j-1)*dx(:,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    tmp = cell(1,nx); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    for i=1:nx ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        tmp(i) = mat2cell(xGridtmp(i,:),1,setup.xPointsDim); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    xGrid = cartprod(tmp); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The grid for numerical itegration ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'tmp_e        = cell(1,model.ne); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'tmp_w        = cell(1,model.ne); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'for i=1:model.ne ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    tmp_e(i) = mat2cell(setup.GH_e,length(setup.GH_e),1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    tmp_w(i) = mat2cell(setup.GH_w,length(setup.GH_w),1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'eGrid = cartprod(tmp_e); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'wGrid = prod(cartprod(tmp_w),2); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Building the residuals ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'f        = zeros(ny+nx,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'ResEuler = zeros(size(xGrid,1),ny+nx); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'yGrid    = zeros(size(xGrid,1),ny); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'for k = 1:size(xGrid,1) ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    x_cu   = transpose(xGrid(k,:)); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    % Getting the states at time t ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = ['    ',char(x(i)),' = model.h0(',num2str(i),',1) + x_cu(', num2str(i),',1);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    % Getting the controls at time t and part of the states at time t+1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    if extendedPerOn == 1';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        [yVec_cu,xVec_cup,diagOut] = policyFunctionEPer(x_cu+model.h0,model,setup.setupEPer);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    else';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        [yVec_cu,xVec_cup]=policy(x_cu+model.h0,model.g0,model.h0,model.derivs,orderApp);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    yGrid(k,:) = transpose(yVec_cu-model.g0);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    % The controls ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
     text = ['    ',char(y(i)),' = yVec_cu(', num2str(i),',1);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    Euler = zeros(ny+nx,1);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    for i=1:length(wGrid)';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        x_cup = xVec_cup-model.h0; ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        % Adding shocks to the states';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        x_cup = x_cup + model.eta*transpose(eGrid(i,:));';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        % Getting controls at time t+1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        if extendedPerOn == 1';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '            yVec_cup = policyFunctionEPer(x_cup+model.h0,model,setup.setupEPer);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        else';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '            yVec_cup = policy(x_cup+model.h0,model.g0,model.h0,model.derivs,orderApp);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        end';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(yp)
     text = ['        ',char(yp(i)),' = yVec_cup(', num2str(i),',1);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '        % Getting the states at time t+1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(xp)
     text = ['        ',char(xp(i)),' = model.h0(',num2str(i),',1) + x_cup(', num2str(i),',1);'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = '        %% Model Specific ';
dlmwrite(saveFile,text,'-append','delimiter','');
DispSymMatrixMatlab_file(strucOfArrays.('f'),'f',saveFile);
text = '        Euler(:,1) = Euler(:,1) + wGrid(i,1)*f(:,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    ResEuler(k,:) = transpose(Euler(:,1));';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');

end